Average Bubble Profile¶
In [1]:
import numpy as np
import math
import matplotlib
#matplotlib.rcParams.update({'font.size': 12})
from matplotlib import gridspec
import matplotlib.pyplot as plt
from collections import deque
import scipy as scp
import scipy.optimize as sco
import scipy.signal as scs
import scipy.special as ssp
import scipy.integrate as sci
import scipy.interpolate as intp
import statistics as stat
import random
from functools import partial
from operator import eq
from itertools import zip_longest, compress, count, islice, groupby, cycle
from labellines import labelLine, labelLines
import os
np_load_old = np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)
Simulation Parameters¶
In [2]:
nLat = 4096
nSims = 390
minSim = 0
simstep = 1
nu = 2.*10**(-3)
lamb = 6; print('lamb = ', lamb)
m2eff = 4. * nu * (- 1. + lamb**2); print('m2eff = ', m2eff)
lenLat = 2 * 50. / np.sqrt(2. * nu); print('lenLat = ', lenLat)
phi0 = 2. * np.pi / 7.; print('phi0 = ', phi0)
alpha = 8.
nCols = 1
phi_initial = np.pi
dt_phi_initial = 0.
mask = 4*phi_initial
############################################################
nyq = nLat//2+1
spec = nyq
hLat = nLat//2
dx = lenLat/nLat; print(dx)
dk = 2.*np.pi/lenLat
dt = dx/alpha
dtout = dt*alpha
light_cone = dtout/dx
In [3]:
titles = [r'$\phi(x)$', r'$\partial_t \phi(x)$', r'$|\nabla \phi(x)|^2$', r'$V(\phi(x))$']
plots_file = '/home/dpirvu/big_plot_file/thin_wall_average_bubble/'
pickle_location = '/home/dpirvu/pickle_location/thin_wall_average_bubble/'
cutout_location = '/gpfs/dpirvu/thin_wall_average_bubble/'
suffix = '_for_phi0{:.4f}'.format(phi0)+'_len{:.4f}'.format(lenLat)+'_lamb{:.4f}'.format(lamb)+'_x'+str(nLat)
def bubbles_file(min, max):
return cutout_location+'bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def centred_bubbles_file(min, max):
return pickle_location+'centred_bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def bubbles_at_rest_file(sim):
return cutout_location+'restbubble_sim'+str(sim)+suffix+'.npy'
def average_bubble(min, max):
return pickle_location+'average_bubble_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def average_of_N_bubbles(Nbubbles):
return pickle_location+'average_of_'+str(Nbubbles)+'_bubbles'+suffix+'.npy'
In [4]:
instanton_location = '/home/dpirvu/inst/instantons/dev/thin_wall_instanton_sim.dat' #/np.sqrt(m2eff)/dx
a = np.genfromtxt(instanton_location)
coleman_profile = np.pi+a[:,0]
xoffset = nLat//8
temp = coleman_profile[len(coleman_profile)//2-xoffset:len(coleman_profile)//2+xoffset]
peaks, _ = scs.find_peaks(temp)
fwhm, height, left_ips, right_ips = scs.peak_widths(temp, peaks, rel_height=0.5)
radius_Coleman_bubble = max(fwhm)/2.
filter_size = radius_Coleman_bubble*dx
print('filter_size = ', filter_size)
plt.plot(np.arange(len(temp))-xoffset, temp)
plt.plot(peaks-xoffset, temp[peaks], "x")
[plt.hlines(height[i], left_ips[i]-xoffset, right_ips[i]-xoffset, color="C2") for i in range(len(fwhm))]
plt.xlabel(r'$r_{\rm E}$'); plt.ylabel(r'$\phi$')
Out[4]:
In [5]:
def V(phi):
return ( -np.cos(phi) + 0.5 * lamb**2. * np.sin(phi)**2. ) * 4. * nu
def dV(phi):
return ( np.sin(phi) + 0.5 * lamb**2. * np.sin(2.*phi) ) * 4. * nu
right_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[np.pi, 2*np.pi], method='bounded')
left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[0, np.pi], method='bounded')
right_left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[2*np.pi, 3*np.pi], method='bounded')
def F(x):
return V(x) - V(phi_initial)
phi_upper_bound = sco.fsolve(F, 5.5)[0]
phi_lower_bound = sco.fsolve(F, 0.5)[0]
phi_lower_bound2 = sco.fsolve(F, -0.5)[0]
print(phi_lower_bound2)
phi_upper_lower_bound = sco.fsolve(F, 6.5)[0]
############################################################
wall_tension, err = sci.quad(lambda x: np.sqrt(2*(V(x) - V(2*np.pi))), np.pi, 2*np.pi); print('wall_tension = ', wall_tension)
epsilon = V(np.pi) - V(2*np.pi); print('epsilon = ', epsilon)
R_coleman = wall_tension/epsilon; print('R_coleman = ', R_coleman)
wall_thkn = right_phi_at_V_max.x / np.sqrt( V(right_phi_at_V_max.x) - V(np.pi) ); print('wall_thkn = ', wall_thkn)
############################################################
matplotlib.rcParams.update({'font.size': 12})
fig = plt.figure(figsize=(7, 4))
plt.plot([i for i in np.arange(0, mask, phi_initial/100)], [V(i) for i in np.arange(0, mask, phi_initial/100)])
plt.plot(phi_upper_bound, V(phi_upper_bound), 'ro')
plt.plot(phi_lower_bound, V(phi_lower_bound), 'co')
plt.plot(phi_lower_bound2, V(phi_lower_bound2), 'co')
plt.plot(phi_upper_lower_bound, V(phi_upper_lower_bound), 'yo')
plt.plot(right_left_phi_at_V_max.x, V(right_left_phi_at_V_max.x), 'ko')
[plt.plot(i*phi_initial, V(i*phi_initial), 'bo') for i in range(5)]
[plt.axvline(left_phi_at_V_max.x + i*phi_initial, color='y', ls='--') for i in range(0,3,2)]
[plt.axvline(right_phi_at_V_max.x + i*phi_initial, color='orange', ls='--') for i in range(0,3,2)]
plt.xlabel(r'$\phi(x)$'); plt.ylabel(r'$V(\phi(x))$'); plt.show()
VISUALISING SIMULATION DATA¶
In [6]:
def add_mask(field_slice, threshold):
return field_slice * [0 if np.abs(i) > threshold else 1 for i in field_slice]
def smoothen(field_slice, sigma, bool):
N = len(field_slice)
if bool:
pbc = [j-1 if j < N/2+1 else N-j+1 for j in range(1, N+1)]
else:
pbc = [j-1 for j in range(1, N+1)]
window = [np.exp(- 0.5 * (x*dx/sigma)**2) / np.sqrt(2*np.pi) / (sigma/dx) for x in pbc]
window = window / sum(window)
spectral_filter = np.fft.fft(window, len(field_slice))
fft_field_slice = np.fft.fft(field_slice, len(field_slice))
smooth_fft_field_slice = [spectral_filter[k] * fft_field_slice[k] for k in range(len(fft_field_slice))]
smooth_field_slice = np.fft.ifft(smooth_fft_field_slice, len(field_slice))
return np.asarray([i.real for i in smooth_field_slice])
def find_slice_peaks(field_slice, peak_threshold):
""" Finds x coordinate of peaks in masked field with mask applied at threshold. """
peak_coord = signal.find_peaks(field_slice, height = peak_threshold)[0].tolist()
if field_slice[-1] > peak_threshold and field_slice[0] > peak_threshold and field_slice[-1] != field_slice[0]:
if field_slice[0] > field_slice[-1] and field_slice[0] > field_slice[1]:
peak_coord.append(0)
elif field_slice[0] < field_slice[-1] and field_slice[-1] > field_slice[-2]:
peak_coord.append(len(field_slice)-1) # this minds potential boundary discontinuities
peak_heights = [field_slice[coord] for coord in peak_coord]
return peak_coord, field_slice.tolist().index(np.max(peak_heights))
def truncateNum(num, decimal_places):
StrNum = str(num)
p = StrNum.find(".") + 1 + decimal_places
return float(StrNum[0:p])
def time_at_fraction(bubble, frac, limit):
T, N = len(bubble), len(bubble[0])
right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
time_list = [t if (right_phi_x[t] <= N*frac) else 0 for t in range(T)]
return next((t for t in time_list[::-1] if t != 0), 0)
def time_at_size(bubble, size, limit):
T, N = len(bubble), len(bubble[0])
right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
time_list = [t if (right_phi_x[t] <= size) else 0 for t in range(T)]
return next((t for t in time_list[::-1] if t != 0), T-1)
def center_bubble(bubble):
bubble0 = bubble[0]
#truncate bubble where it expands relativistically
limit = phi_upper_bound
tdecap = time_at_fraction(bubble0, 0.9, np.floor(limit))
tmax = time_at_fraction(bubble0, 0.95, 4.5)
# tmax = time_at_fr(bubble0, 0.01, limit)
T, N = len(bubble0), len(bubble0[0])
#rotate bubble such that the relativistic bit is perfectly centered
fld = smoothen(bubble0[tdecap], 1, True)
kafslj = next(iter([x for x in range(N-1) if fld[x] > limit]))
if kafslj < 30:
fld = smoothen(bubble0[tdecap], dx, True)
vals = [x if (fld[x] > limit and fld[x+1] > limit) else 0 for x in range(N-1)]
bubbles = [list(g) for k, g in groupby(vals, lambda x: x != 0) if k]
maxi, imax, keep = 0, 0, 0
for i in range(len(bubbles)-1):
maxi = bubbles[i+1][0]-bubbles[i][-1]
if maxi > imax:
imax = maxi
keep = i
# first_zero = next((i for i, x in enumerate(vals) if x == 0.), 0) # where bubble begins
first_zero = int(np.mean([bubbles[keep][0], bubbles[keep][-1]]))
fld = np.roll(bubble0[tdecap], -first_zero)
target_peak = int(round(np.mean([x for x in range(N-1) if (fld[x] > limit and fld[x+1] > limit)]))) # center of bubble
# print(first_zero, target_peak)
angle = int(N//2) - (target_peak + first_zero) # rotation angle needed
bubble = np.asarray([[np.roll(bubble[col][t], angle) for t in range(tmax)] for col in range(len(bubble))])
else:
bubble = np.asarray([bubble[col][:tmax] for col in range(len(bubble))])
# check which sides the COM travels
bubble0 = bubble[0]
tcheck = time_at_fraction(bubble0, 0.02, 2*phi_initial)
# print(tcheck)
fld = [smoothen(bubble0[t], dx, True) for t in range(tcheck-100, tcheck, 2)]
# tv = [np.nanmean([x for x in range(N-1) if (fld[t][x] > 2*phi_initial and fld[t][x+1] > 2*phi_initial)]) for t in range(len(fld))]
tv = [np.nanmean([x for x in range(N) if fld[t][x] > 2*phi_initial]) for t in range(len(fld))]
decr = np.sum([1 for i, j in zip(tv, tv[1::]) if i > j])
incr = np.sum([1 for i, j in zip(tv, tv[1::]) if i < j])
# print(incr, decr)
return bubble, decr - incr #np.nanmean(tv) - N//2
def multiply_bubble(bubble, dir, fold):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
bubble = [np.tile(bubble[col], fold) for col in range(len(bubble))] # multiplies bubbles so tail is kept without pbc
TT, NN = len(bubble[0]), len(bubble[0][0])
for t in range(TT):
a, b = int((TT-t)*light_cone) + N - 200, int((TT-t)*light_cone/2) - N//4
x1, x2 = np.arange(a, NN), np.arange(b)
if dir < 0:
x1, x2 = x1 - a, x2 - (b-NN)
for x in np.append(x1, x2):
if 0 <= x < NN:
bubble[0][t][x] = phi_initial
return np.asarray(bubble)
In [7]:
#all_data1, sims_to_keep1 = np.load(bubbles_file(minSim, 120))
#all_data2, sims_to_keep2 = np.load(bubbles_file(120, 240))
#all_data3, sims_to_keep3 = np.load(bubbles_file(240, 340))
#all_data4, sims_to_keep4 = np.load(bubbles_file(340, nSims))
#all_data = np.concatenate((all_data1, all_data2, all_data3, all_data4)).tolist()
#sims_to_keep = np.concatenate((sims_to_keep1, sims_to_keep2, sims_to_keep3, sims_to_keep4)).tolist()
#del all_data1, all_data2, all_data3, all_data4, sims_to_keep1, sims_to_keep2, sims_to_keep3, sims_to_keep4
# if trying to save the kernel:
all_data, sims_to_keep = np.load(bubbles_file(340, nSims))
In [8]:
normal = [phi_initial, np.mean([sim[1][0] for sim in all_data]), np.mean([sim[2][0] for sim in all_data]), V(phi_initial), 0.5*np.mean([sim[1][0] for sim in all_data])**2]
print(normal)
In [9]:
#bubble_data = []
#for sim in range(len(all_data))[1:1:]:
# print(sim)
# col = 0
# simulation = all_data[sim]
# centered, dir = center_bubble(simulation)
# mult = multiply_bubble(centered, dir, 2)
# fig, ax = plt.subplots(1, 3, figsize = (12, 3))
# im0 = ax[0].imshow(simulation[col], aspect='auto', interpolation='none', origin='lower')
# im1 = ax[1].imshow(centered[col], aspect='auto', interpolation='none', origin='lower')
# im2 = ax[2].imshow(mult[col], aspect='auto', interpolation='none', origin='lower')
# clb = plt.colorbar(im0, ax = ax[0]); clb = plt.colorbar(im1, ax = ax[1]); clb = plt.colorbar(im2, ax = ax[2]); plt.show()
# bubble_data.append(mult)
#np.save(centred_bubbles_file(minSim, nSims), bubble_data)
#check_multiplication = np.load(centred_bubbles_file(minSim, nSims))
#for pic in check_multiplication:
# plt.figure(0); plt.imshow(pic, aspect='auto', interpolation='none', origin='lower'); plt.show()
In [10]:
def rapidity(vel):
return np.arctanh(vel)
def gamma(v):
return ( 1. - v**2. )**(-0.5)
def get_bubble_radius(bubble_slice, filter):
bubble_slice = smoothen(bubble_slice, filter, True)
results_half = signal.peak_widths(bubble_slice, signal.find_peaks(bubble_slice)[0], rel_height=0.5)
return max(results_half[0])/2.
def tanh_pos(x, r0L, r0R, dr, vL, vR):
wL, wR = dr/gamma(vL), dr/gamma(vR)
return ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi
def tanh_neg(x, r0L, r0R, dr, vL, vR):
wL, wR = dr/gamma(vL), dr/gamma(vR)
return - ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi
def tanh_fit(bubble_slice, axis, prior):
# plt.plot(axis, bubble_slice, 'ro', axis, [tanh(r, *best_tanh) for r in axis], 'g'); plt.show()
bounds = ((axis[0], 0, 0, 0, 0), (0, axis[-1], axis[-1], 1, 1))
if prior is not None:
if np.mean(bubble_slice) > phi_initial:
return sco.curve_fit(tanh_pos, axis, bubble_slice, p0=prior, bounds=bounds)[0]
else:
return sco.curve_fit(tanh_neg, axis, bubble_slice, p0=prior, bounds=bounds)[0]
else:
if np.mean(bubble_slice) > phi_initial:
return sco.curve_fit(tanh_pos, axis, bubble_slice, bounds=bounds)[0]
else:
return sco.curve_fit(tanh_neg, axis, bubble_slice, bounds=bounds)[0]
def hyperbola1(t, a, b, c):
return np.sqrt(c + (t - b)**2) + a
def hyperbola2(t, d, e, f):
return - np.sqrt(f + (t - e)**2) + d
def hypfit_test(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
try:
if rr[0] <= rr[-1]:
fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
else:
fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
return True
except (RuntimeError, TypeError):
return False
def hypfit(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
if rr[0] <= rr[-1]:
fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
traj1, traj2 = hyperbola1(tt, *fit1), hyperbola2(tt, *fit2)
else:
fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
traj1, traj2 = hyperbola2(tt, *fit1), hyperbola1(tt, *fit2)
# print('hyperbolic trajectories fit: ', fit1, fit2)
# plt.plot(rr, tt, 'g-', qq, tt, 'r-', traj1, tt, 'y:', traj2, tt, 'b:') # plot the equation using the fitted parameters
return tt, traj1, traj2
def get_velocities(tt, rrwallfit, llwallfit):
#uu = wall travelling along with the com; vv = wall travelling against the com; aa = centre of mass velocity; bb = instantaneous velocity of wall
uu = intp.splev(tt, intp.splrep(tt, rrwallfit), der=1)
vv = intp.splev(tt, intp.splrep(tt, llwallfit), der=1)
for i in range(len(uu)):
if np.abs(uu[i]) > 1:
uu[i] = np.sign(uu[i])*(1-np.abs(np.abs(uu[i])-1))
if np.abs(vv[i]) > 1:
vv[i] = np.sign(vv[i])*(1-np.abs(np.abs(vv[i])-1))
# print(uu)
# print(vv)
for i in range(len(uu)):
if str(uu[i]) == 'nan' :
uu[i] = - np.sign(vv[-1])*0.999
if str(vv[i]) == 'nan' :
vv[i] = - np.sign(uu[-1])*0.999
# print(uu)
# print(vv)
# if all(np.abs(i) > 1 for i in uu):
# uu = np.asarray([np.sign(i)*(1-np.abs(np.abs(i)-1)) for i in uu])
# if all(np.abs(i) > 1 for i in vv):
# vv = np.asarray([np.sign(i)*(1-np.abs(np.abs(i)-1)) for i in vv])
aa = ( 1 + uu*vv - np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / ( uu + vv)
bb = (-1 + uu*vv + np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / (-uu + vv)
return uu, vv, aa, bb#, ch1, ch2
def velocity(bubble, bool):
T, N = len(bubble), len(bubble[0])
limit = phi_upper_bound - 0.5 #phi_upper_bound# - 0.5*np.abs(phi_upper_bound-np.floor(phi_upper_bound))
prior = None
data_list = []
err = 0
# find where fit can begin
window = int(N//10)
if window > 100: window = 100
tf = T-1
endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
if len(endSlice) == 0:
return 'next'
while endSlice[0] - window < 0 or endSlice[-1] + window >= N:
tf -= 2
endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
if len(endSlice) == 0:
return 'next'
for t in range(tf, 0, -1):
if len(data_list) == 0:
peaks = endSlice
target = int(np.round(np.nanmean(peaks)))
coord_list = np.arange(peaks[0] - window, peaks[-1] + window)
else:
xrange = np.arange(data_list[-1][0]-window, data_list[-1][0]+window+1)
if all(0 <= i < N for i in xrange) and any(np.cos(bubble[t][i]) > np.cos(limit) for i in xrange):
peaks = [i for i in xrange if np.cos(bubble[t][i]) > np.cos(limit)]
else:
break
target = int(np.round(np.nanmean(peaks)))
coord_list = np.arange(target - int(np.round(np.abs(data_list[-1][1]))) - window, target + int(np.round(np.abs(data_list[-1][2]))) + window)
coords = np.asarray([(bubble[t][i%N], int(i-target)) for i in coord_list])
if err < 5:
try:
r0L, r0R, dr, vL, vR = tanh_fit(coords[:,0], coords[:,1], prior)
prior = r0L, r0R, dr, vL, vR
data_list.append([target, r0L, r0R, int(t)])
except (RuntimeError, ValueError, TypeError):
prior = None
data_list.append([target, r0L + np.sign(data_list[-1][1]-data_list[0][1])*light_cone*(T/N), r0R + np.sign(data_list[-1][2]-data_list[0][2])*light_cone*(T/N), int(t)])
err += 1
else:
break
data_list = np.asarray(data_list[::-1])
targets, r0Ls, r0Rs, time_list = data_list[:,0], data_list[:,1], data_list[:,2], data_list[:,-1]
# print('tmin, tmax', time_list[0], time_list[-1])
# get direction of bubble
radius_diff = np.mean([np.abs(data_list[i,1]) - np.abs(data_list[i,2]) for i in range(len(data_list)//4)])
if radius_diff > 0: # if average difference in radius is positive then left radius is on average higher than right radius i.e. bubble travels to the right
rr, ll = targets + r0Rs, targets + r0Ls
else:
rr, ll = targets + r0Ls, targets + r0Rs
if bool:
fig, ax0 = plt.subplots(1, 1, figsize = (5, 4))
ax0.plot(rr, time_list, 'b', ll, time_list, 'y', linewidth='3')
ax0.set_xlabel('x'); ax0.set_ylabel('t')
# get velocities from derivative of best fit to wall trajectory
trunc = 0
# save copies
time_list_copy, rr_copy, ll_copy = time_list, rr, ll
# as lons as it needs..
while(True):
# try to fit walls
if hypfit_test(time_list_copy, rr_copy, ll_copy):
time_list, rrwallfit, llwallfit = hypfit(time_list_copy, rr_copy, ll_copy)
uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)
break
# if fit fails, shorten walls
else:
trunc += 25
if len(time_list_copy) > 100 + trunc*2:
time_list_copy, rr_copy, ll_copy = time_list_copy[trunc:-trunc:], rr_copy[trunc:-trunc:], ll_copy[trunc:-trunc:]
if trunc > 200 or len(time_list_copy) < 100:
return 'next'
continue
# if hypfit_test(time_list, rr, ll):
# time_list, rrwallfit, llwallfit = hypfit(time_list, rr, ll)
# else:
# return 'next'
# uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)
# if bool:
# fig, [ax0, ax1] = plt.subplots(1, 2, figsize = (15, 4))
# ax0.plot(rr[-len(time_list):], time_list, 'b', ll[-len(time_list):], time_list, 'y', rrwallfit, time_list, 'r:', llwallfit, time_list, 'r:', linewidth='3')
# ax0.set_xlabel('x'); ax0.set_ylabel('t')
# ax1.plot(time_list, uu, 'b:', time_list, vv, 'y:')
# ax1.plot(time_list, aa, 'r', label=r'v COM')
# ax1.plot(time_list, bb, 'g', label=r'v walls')
# ax1.axhline(0, color='darkgray', ls=':')
# ax1.set_xlabel('t'); ax1.set_ylabel('v(t)'); ax1.legend(); plt.show()
# plt.savefig(plots_file+'velocity_profile'+str(random.randrange(100))+suffix+'.png');
list = np.abs(aa - bb)
if any(str(i) != 'nan' for i in list):
return -aa[list.tolist().index(np.nanmin(list))]
else:
return 'next'
In [11]:
def fold(beta):
if np.abs(beta) > 0.8:
return 4
else:
return 3
def add_velocities(v1,v2):
return (v1 + v2) / (1. + v1*v2)
def cut_out_bubble(bubble):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
ext = int(N//8)
if ext < 100: ext = 200
xmin = next((i for i in range(N) if bubble0[-1][i] > phi_upper_bound), 0) - ext
xmax = next((i for i in range(N)[::-1] if bubble0[-1][i] > phi_upper_bound), N-1) + ext
if xmin < 0: xmin = 0
if xmax >= N: xmax = N-1
t, tmin = 0, 0
while not any(np.cos(bubble0[t][i]) > np.cos(3.17) for i in range(xmin, xmax+1)):
tmin = t
t += 1
return [[[bubble[col][t][x] for x in range(xmin, xmax+1)] for t in range(tmin, T)] for col in range(len(bubble))]
def center_bubble_again(bubble, kind):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
ext = int(N//8)
if ext < 100: ext = 200
if kind < 0:
x1 = 0
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
while np.cos(np.nanmean(peaks1)) < np.cos(3.19) and x1 < N-1:
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
x1 += 1
if len(peaks2) > 0:
tmax = peaks2[-1]
else:
return 'next'
# print('tmax, x1', tmax, x1)
xmin = x1 - ext
if xmin < 0: xmin = 0
xmax = next((i for i in range(N-1, xmin, -1) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), N-1)
else:
x1 = N-1
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
while np.cos(np.nanmean(peaks1)) < np.cos(3.19) and x1 > 0:
slice = [bubble0[i][x1] for i in range(T)]
peaks1 = [i for i in slice if i < mask]
peaks2 = [i for i in range(T) if slice[i] < mask]
x1 -= 1
if len(peaks2) > 0:
tmax = peaks2[-1]
else:
return 'next'
# print('tmax, x1', tmax, x1)
xmax = x1 + ext
if xmax > N-1: xmax = N-1
xmin = next((i for i in range(xmax) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), 0)
t, tmin = 0, 0
while not any([(mask > bubble0[t][i] and np.cos(bubble0[t][i]) > np.cos(phi_initial)) for i in range(xmin, xmax+1)]):
tmin = t
t += 1
return [[[bubble[col][t][x] if bubble[col][t][x] != mask else normal[col] for x in range(xmin, xmax+1)] for t in range(tmin, tmax+1)] for col in range(len(bubble))]
def boost(vCOM, ga, x, t):
return ga * (x - vCOM*t)
def interpolate_bubble(bubble, res):
NT, N = len(bubble), len(bubble[0])
bubble = np.asarray(bubble)
t, x = np.arange(NT), np.arange(N)
tt, xx = np.meshgrid(t, x)
f = intp.interp2d(t, x, bubble[tt, xx], kind = 'quintic')
T, X = np.arange(0, NT, 1/res), np.arange(0, N, 1/res)
return f(T, X).T
def boost_bubble(bubble, vCOM, res):
bubble0 = bubble[0]
T, N = len(bubble0), len(bubble0[0])
limit = np.ceil(phi_upper_bound)
ga = gamma(vCOM)
t0 = int(time_at_size(bubble0, 100, limit)) # arbitrary choice of bubble centre
x0 = np.mean([i for i in range(N) if mask > bubble0[t0][i] > limit])
print(t0, x0)
T0, X0 = boost(vCOM, ga, t0, x0), boost(vCOM, ga, x0, t0)
deltaT, deltaX = np.abs(T0-t0), np.abs(X0-x0)
bubble = [interpolate_bubble(bubble[col], res) for col in range(len(bubble))]
new_bubble = [[[mask for x in range(N)] for t in range(T)] for col in range(len(bubble))]
for t in range(T):
for x in range(N):
tt = int((boost(vCOM, ga, t, x) + np.sign(T//2-t0) * deltaT)*res)
xx = int((boost(vCOM, ga, x, t) + np.sign(N//2-x0) * deltaX)*res)
if (T*res > tt >= 0 and N*res > xx >= 0):
for col in range(len(bubble)):
new_bubble[col][t][x] = bubble[col][tt][xx]
return new_bubble, x0 - N//2
In [12]:
def the_whole_shebang():
resolution = 3
for sim in range(len(all_data))[10::]:
exists = os.path.exists(bubbles_at_rest_file(sim))
if not exists:
print(sim)
bubble, dir = center_bubble(all_data[sim])
beta = velocity(bubble[0], True)
if beta == 'next':
print('sim '+str(sim)+' skipped')
continue
elif np.abs(beta) > 0.9:
beta = np.sign(beta)*0.9
bubble = multiply_bubble(bubble, dir, fold(beta))
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.show()
res, vtotal = 0, 0
bool1, bool2 = True, True
while np.abs(beta) > 0.05:
bool2 = True
bubble, dir = boost_bubble(bubble, beta, resolution)
bubble = center_bubble_again(bubble, dir)
if bubble == 'next':
print('sim '+str(sim)+' interrupted')
bool1 = False
break
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.show()
prebeta = beta
vtotal = add_velocities(vtotal, beta)
print('step done ', beta)
min_bubble = cut_out_bubble(bubble)
fig = plt.figure(figsize=(5,5))
plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
plt.show()
beta = velocity(min_bubble[0], True)
if beta == 'next':
print('sim '+str(sim)+' interrupted')
if np.abs(prebeta) >= 0.8:
bool1 = False
break
elif np.abs(beta) > 0.9:
beta = np.sign(beta)*0.9
# res += 1
if not bool2:
if np.abs(beta) > 0:
min_bubble = cut_out_bubble(bubble)
vtotal = beta
bool2 = True
if bool1 and bool2:
print('sim', sim, 'total vel ', vtotal)
np.save(bubbles_at_rest_file(sim), [min_bubble, vtotal])
return
In [13]:
#the_whole_shebang()
Average bubble¶
In [14]:
def cut_out_coordinates(bubble0):
T, N = len(bubble0), len(bubble0[0])
limit = np.floor(phi_upper_bound)
size = 150 # do not change
tmin = T-1
for t in range(T):
slice = smoothen(bubble0[t], dx, False)
right_phi = [x for x in range(N-1) if (slice[x] >= limit and slice[x+1] >= limit)]
if len(right_phi) >= size:
tmin = t
break
if tmin == T-1:
return None
if stat.stdev(right_phi) > 100:
return None
xcen = int(np.round(np.nanmean(right_phi)))
xcen_list = [np.round(np.nanmean([x for x in range(N) if bubble0[t][x] >= limit])) for t in range(tmin-50,tmin,2)]
if stat.stdev(xcen_list) > 50:
return None
# fig = plt.figure(figsize=(4,4))
# plt.plot(xcen, tmin, 'ro')
# plt.imshow(bubble, aspect='auto', interpolation='none', origin='lower')
# plt.show()
ex1 = np.arange(tmin)[::-1]
ex2 = np.arange(tmin, T)
ex3 = np.arange(xcen)[::-1]
ex4 = np.arange(xcen, N)
if len(ex1) > 500: ex1 = ex1[:500]
if len(ex2) > 500: ex2 = ex2[:500]
if len(ex3) > 500: ex3 = ex3[:500]
if len(ex4) > 500: ex4 = ex4[:500]
return ex1, ex2, ex3, ex4
def stack_bubbles(bubble_data):
extent_list, bad_bubbles = [], []
for sim in range(len(bubble_data)):
coords = cut_out_coordinates(bubble_data[sim][0])
if coords is not None:
extent_list.append(coords)
else:
bad_bubbles.append(sim)
extent_list = np.asarray(extent_list)
print('bad bubbles', bad_bubbles)
bubble_data = [bubble_data[i] for i in range(len(bubble_data)) if i not in bad_bubbles]
nCols = 5
for sim in range(len(bubble_data)):
bubble_data[sim].append([0.5*np.asarray(i)**2 for i in bubble_data[sim][1]])
uper_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
upper_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
lower_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
lower_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
av_mat, av_err_mat = [], []
for col in range(nCols):
average_matrix, average_error_matrix = [[],[],[],[]], [[],[],[],[]]
coord_lists = [uper_right_bubble_list[col], upper_left_bubble_list[col], lower_right_bubble_list[col], lower_left_bubble_list[col]]
for i in range(4):
max_lines, max_cols = 500, 500
for line in range(max_lines):
average_matrix[i].append([normal[col] for n in range(max_cols)])
average_error_matrix[i].append([0. for n in range(max_cols)])
for num_line in range(max_lines):
for num_col in range(max_cols):
meas = []
for simulation in coord_lists[i]:
if len(simulation) > num_line and len(simulation[0]) > num_col:
meas.append(simulation[num_line][num_col])
average_matrix[i][num_line][num_col] = np.mean(meas)
if len(meas) != 0.:
average_error_matrix[i][num_line][num_col] = stat.stdev(meas)/len(meas)
av_mat.append(average_matrix)
av_err_mat.append(average_error_matrix)
return [av_mat, av_err_mat], bad_bubbles
def restore_average_bubble(stacked_bubble_data):
whole_bubbles = []
for av_bub in stacked_bubble_data:
whole_bubbles.append([])
for col in range(len(av_bub)):
top = np.concatenate((np.flip(np.asarray(np.flip(av_bub[col][0],1)).transpose(),1), np.flip(np.asarray(av_bub[col][1]).transpose(),0)), axis=1)
bottom = np.concatenate((np.flip(np.asarray(av_bub[col][2]).transpose(),1), np.flip(np.asarray(np.flip(av_bub[col][3],0)).transpose(),1)), axis=1)
whole_bubbles[-1].append(np.concatenate((top, bottom), axis=0))
return np.asarray(whole_bubbles)
In [15]:
bubble_list = []
for sim in range(200):#len(all_data)):
bool = os.path.exists(bubbles_at_rest_file(sim))
if bool:
bubble_list.append(np.load(bubbles_at_rest_file(sim))) #deboosted_bubbles2 #deboosted_bubbles_loop_all_sims7
bubble_list = np.asarray(bubble_list)
bubble_list_save = bubble_list[:,0]
bubble_velocities = bubble_list[:,1]
In [16]:
#if testing
#normal = [phi_initial, 0, 0, V(phi_initial), 0]
#either compute average
stacked_bubble_data, bad_bubbles = stack_bubbles(bubble_list_save)
ab, error_ab = restore_average_bubble(stacked_bubble_data)
np.save(average_bubble(minSim, nSims), [ab, error_ab])
#or load it
#ab, error_ab = np.load(average_bubble(minSim, nSims))
#separate columns
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
error_bubble, error_mom_bubble, error_ge_bubble, error_pe_bubble, error_ke_bubble = error_ab
In [17]:
#plot discarded bubbles
for i in range(len(bubble_list_save)):
if i in bad_bubbles:
fig, ax = plt.subplots(1, 1, figsize=(7,5))
im = ax.imshow(bubble_list_save[i][0], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax)
plt.show()
In [18]:
bubble_list_save = [bubble_list_save[i] for i in range(len(bubble_list_save)) if i not in bad_bubbles]
bubble_velocities = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i not in bad_bubbles]
check = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i in bad_bubbles]
print('Total number of bubbles available: ', len(bubble_list_save))
In [19]:
for i in bubble_list_save:
fig, ax = plt.subplots(1, len(i), figsize=(25,5))
for col in range(len(i)):
im = ax[col].imshow(i[col], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax[col])
plt.show()
In [20]:
#plot average bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble
fig, ax = plt.subplots(1, len(ab), figsize=(25, 5))
for col in range(len(ab)):
im = ax[col].imshow(ab[col], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax[col])
plt.savefig(plots_file+'average_bubble'+suffix+'.png')
#plot distribution of boosts
fig, ax = plt.subplots(1, 3, figsize=(15,5))
ax[0].hist(bubble_velocities, bins=len(bubble_velocities))
ax[1].hist([gamma(i) for i in bubble_velocities], bins=len(bubble_velocities))
ax[2].hist([rapidity(i) for i in bubble_velocities], bins=len(bubble_velocities))
[ax[i].set_ylabel('Count') for i in range(len(ax))]
ax[0].set_xlabel(r'$v$'); ax[1].set_xlabel(r'$\gamma$'); ax[2].set_xlabel(r'$w$');
plt.savefig(plots_file+'boost_distrib'+suffix+'.png');
In [21]:
def coleman_match(bubble, colemanSoln):
minSum = np.sum(colemanSoln)**2.
tColeman = 0
for t in range(len(bubble)):
slice = smoothen(bubble[t], 2*dx, True)
a = sum([(slice[x] - colemanSoln[x])**2. for x in range(len(slice))])
if a < minSum:
minSum, tColeman = a, t
return tColeman
def plot_coleman_match(bubble, error_bubble, colemanSoln):
T, N = len(bubble), len(bubble[0])
colemanSoln = colemanSoln[len(colemanSoln)//2-N//2:len(colemanSoln)//2+N//2]
tcoleman = coleman_match(bubble, colemanSoln)
slice, err_slice = bubble[tcoleman], error_bubble[tcoleman]
fig = plt.figure(figsize=(15,5))
plt.errorbar(np.arange(N), slice, yerr = err_slice, color = 'k', marker='o', ms=1, alpha = 0.8, ecolor = 'green', label = 'Average bubble')
plt.plot(np.arange(N), colemanSoln, 'r-', label = 'Best match t = '+str(tcoleman)); plt.legend()
plt.savefig(plots_file+'average_bubble_profile'+suffix+'.png')
return tcoleman
def bubble_total_energy(KED, GED, PED):
T, N = len(KED), len(KED[0])
KE = [sum(KED[t]) for t in range(T)]
GE = [sum(GED[t]) for t in range(T)]
PE = [sum(PED[t]) for t in range(T)]
TE = [sum(x) for x in zip(KE, GE, PE)]
return KE, GE, PE, TE
def plot_total_energy_density(KED, GED, PED, tcoleman):
energies = [KED, GED, PED, KED + GED + PED]
titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']
fig, ax = plt.subplots(2, 2, figsize=(15,10))
N = len(KED[0])
for i in range(4):
if i in [0,1]: j = 0
else: j = 1
im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
ax[i%2][j].plot(np.arange(N), np.ones(N)*tcoleman, 'r-')
plt.savefig(plots_file+'energy_densities_from_average'+suffix+'.png');
return
def plot_total_energy(KED, GED, PED, tcoleman):
energies = bubble_total_energy(KED, GED, PED)
titles = ['KE', 'GE', 'PE', 'TE']
ls = ['b-', 'g-', 'y-', 'r-']
N = len(energies[0])
fig = plt.figure(figsize=(8,5))
for i in range(len(energies)):
plt.plot(np.arange(N), energies[i], ls[i], label = titles[i])
plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
plt.axvline(tcoleman, color = 'darkgray', ls = '--')
plt.xlabel('$t$'); plt.legend(); plt.savefig(plots_file+'energies'+suffix+'.png'); plt.show()
return
def plot_energy_dens_at_t(KED, GED, PED, timeslice):
TED = [sum(x) for x in zip(KED, GED, PED)]
energies = [KED, GED, PED, TED]
titles = ['KED', 'GED', 'PED', 'TED']
ls = ['b-', 'g-', 'y-', 'r-']
N = len(KED[0])
fig = plt.figure(figsize=(8,8)); plt.subplots_adjust(hspace=.0)
gs = gridspec.GridSpec(len(energies), 1, height_ratios=[1]*len(energies))
ax = [[]]*len(energies); ax[0] = plt.subplot(gs[0])
for i in range(1, len(ax)): ax[i] = plt.subplot(gs[i], sharex = ax[i-1])
for i in range(len(ax)):
ax[i].plot(np.arange(N), energies[i][timeslice], ls[i], label = titles[i])
if i == 0:
ax[i].axhline(np.pi, color = 'darkgray', ls='--')
else:
ax[i].axhline(np.mean(energies[i][timeslice]), color = 'darkgray', ls='--')
ax[i].legend()
ax[-1].set(xlabel='$x$')
plt.savefig(plots_file+'energy_densities_at_t'+str(timeslice)+suffix+'.png'); plt.show()
return
In [22]:
TColeman = plot_coleman_match(bubble, error_bubble, coleman_profile)
plot_total_energy_density(ke_bubble, ge_bubble, pe_bubble, TColeman)
plot_total_energy(ke_bubble, ge_bubble, pe_bubble, TColeman)
#plot_energy_dens_at_t(ke_bubble, ge_bubble, pe_bubble, TColeman)
In [23]:
def GradEnDen(field):
Nx = len(field)
NNx = Nx//2+1
fft_field = scp.fft.fft(field)
for i in range(len(fft_field)):
fft_field[i] = 1j*i*dk * fft_field[i]# / Nx
real_field = scp.fft.ifft(fft_field)
return [0.5*i.real**2 for i in real_field]
def plot_derived_total_energy_density(field, momentum, filter, t_coleman):
smooth_field = [smoothen(slice, filter, True) for slice in field]
smooth_momentum = [smoothen(slice, filter, True) for slice in momentum]
KED = [0.5*slice**2. for slice in smooth_momentum]
GED = [GradEnDen(slice) for slice in smooth_field]
PED = [V(slice) for slice in smooth_field]
TED = [[KED[t][x]+GED[t][x]+PED[t][x] for x in range(len(KED[0]))] for t in range(len(KED))]
energies = [KED, GED, PED, TED]
titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']
fig, ax = plt.subplots(2, 2, figsize=(15,10))
N = len(smooth_field[0])
for i in range(4):
if i in [0,1]: j = 0
else: j = 1
im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower')
clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
ax[i%2][j].plot(np.arange(N), np.ones(N)*t_coleman, 'r-')
plt.savefig(plots_file+'energy_densities'+suffix+'.png');
return
def plot_TE_with_Nbubbles(tcoleman):
bubble_averages_list = []
for i in range(Nbubs):
bool = os.path.exists(average_of_N_bubbles(i+1))
if bool:
bubble_averages_list.append(np.load(average_of_N_bubbles(i+1)))
fig = plt.figure(figsize=(8,6))
for bubble_average in bubble_averages_list:
ab, error_ab, length = bubble_average
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
energies = bubble_total_energy(ke_bubble, ge_bubble, pe_bubble)
plt.plot(np.arange(len(energies[-1])), energies[-1], label = length)
plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
plt.axvline(tcoleman, color = 'darkgray', ls = '--')
plt.xlabel('$t$'); plt.legend()
plt.savefig(plots_file+'TE_vs_nbubbles1'+suffix+'.png');
fig = plt.figure(figsize=(8,6))
for bubble_average in bubble_averages_list:
ab, error_ab, length = bubble_average
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
GED, PED, KED = ge_bubble[tcoleman], pe_bubble[tcoleman], ke_bubble[tcoleman]
TED = GED + PED + KED
plt.plot(length, sum(TED), 'ro')#, np.mean(TED))
plt.xlabel('$N$'); plt.ylabel('$E$'); plt.legend()
plt.savefig(plots_file+'TE_vs_nbubbles2'+suffix+'.png');
return
def plot_TE_with_filter(tcoleman):
bubble_average = np.load(average_of_N_bubbles(Nbubs))
ab, error_ab, length = bubble_average
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
slice = bubble[tcoleman]
mom_slice = mom_bubble[tcoleman]
fig = plt.figure(figsize=(8,6))
for filter in np.linspace(dx, 1000*dx, 200):
smooth_field = smoothen(slice, filter, True)
smooth_momentum = smoothen(mom_slice, filter, True)
KED = 0.5*smooth_momentum**2.
GED = GradEnDen(smooth_field)
PED = V(smooth_field)
TED = KED + GED + PED
plt.plot(filter, np.sum(TED), 'ro', ms=2)#, 'o', label='{:.3f}'.format(filter))
plt.axhline(0, ls='--', color='darkgray')
plt.xlabel(r'$\sigma$'); plt.ylabel('$E$'); plt.legend()
plt.savefig(plots_file+'TE_vs_filter'+suffix+'.png');
return
In [24]:
plot_derived_total_energy_density(bubble, mom_bubble, 5, TColeman)
In [25]:
#average x bubbles at a time; save
Nbubs = len(bubble_list_save)
step = 4
for i in range(0, Nbubs, step):
qq = len(bubble_list_save)
bubble_puzzle, _ = stack_bubbles(bubble_list_save)
ab, error_ab = restore_average_bubble(bubble_puzzle)
np.save(average_of_N_bubbles(qq), [ab, error_ab, qq])
for i in range(step):
if len(bubble_list_save) > 0:
del bubble_list_save[0]
In [26]:
plot_TE_with_Nbubbles(TColeman)
In [27]:
#Nbubs = 62
plot_TE_with_filter(TColeman)
In [ ]: